EEID-Microbiome-TimeSeries

Packages

Packages

library(tidyverse)
library(ggplot2)
library(ggthemes)
library(RColorBrewer)
library(lme4)
library(qiime2R)
library(ecodist)
library(vegan)
library(ggpubr)

Bacterial 16S

Base Data Read-in

Reading in Data

Metadata

Metadata16S = read_tsv("./16S/Metadata_NSFEEID_16SRuns_1234Merged_plusExpData.txt")

ReadCount16S = read_tsv("16S/WithPlasmid/PlasmidWork_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_merged1234/CorrectedRead_Counts.txt")

OTU Table

#OTU Tables
AdNewt_Rarified_Tableqza = read_qza("./16S/WithPlasmid/CoreMet_CopCorTrueAbundEst_r5000_Phylo/rarefied_table.qza") 

AdNewt_Rarified_Table = as.data.frame(AdNewt_Rarified_Tableqza$data)

AdNewt_UnRarified_Tableqza = read_qza("./16S/WithPlasmid/PlasmidWork_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_merged1234/CopCor_True_abund_estimate_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_TempDose_Exp_merged1234.qza")

AdNewt_UnRarified_Table = as.data.frame(AdNewt_UnRarified_Tableqza$data)

Taxonomy<-read_qza("./16S/WithPlasmid/EEID_16S_1234_Taxonomywplasmid.qza")$data %>% parse_taxonomy()

Alpha Diversity

#Alpha

AdNewt_PD16S = read_qza("./16S/WithPlasmid/CoreMet_CopCorTrueAbundEst_r5000_Phylo/faith_pd_vector.qza")
AdNewt_PD16S <- tibble::rownames_to_column(AdNewt_PD16S$data, "SampleID")

AdNewt_sOTU16S = read_qza("./16S/WithPlasmid/CoreMet_CopCorTrueAbundEst_r5000_Phylo/observed_features_vector.qza")
AdNewt_sOTU16S <- tibble::rownames_to_column(AdNewt_sOTU16S$data, "SampleID")

AdNewt_Evenness16S = read_qza("./16S/WithPlasmid/CoreMet_CopCorTrueAbundEst_r5000_Phylo/evenness_vector.qza")
AdNewt_Evenness16S <- tibble::rownames_to_column(AdNewt_Evenness16S$data, "SampleID")

AdNewt_Shannon16S = read_qza("./16S/WithPlasmid/CoreMet_CopCorTrueAbundEst_r5000_Phylo/shannon_vector.qza")
AdNewt_Shannon16S <- tibble::rownames_to_column(AdNewt_Shannon16S$data, "SampleID")

Antifungal predictions

#Antifungal 
AntiFungal_Bsal = read_tsv("./16S/WithPlasmid/Antifungal_Predictions_EEID_Bsalonly_nonrar/Metadata_Antifungal_Predictions.txt") %>%
  select(SampleID,AntiFungal_Richness,TotalAntiFungal) %>%
  rename(AntiFungal_Richness_Bsal = AntiFungal_Richness) %>%
  rename(TotalAntiFungal_Bsal = TotalAntiFungal)

AntiFungal_Full = read_tsv("./16S/WithPlasmid/Antifungal_Predictions_EEID_Full_nonrar/Metadata_Antifungal_Predictions.txt") %>%
  select(SampleID,AntiFungal_Richness,TotalAntiFungal)

Beta Diversity

#Beta diversity

AdNewt_Jac16Smx<-read_qza("./16S/WithPlasmid/CoreMet_CopCorTrueAbundEst_r5000_Phylo/jaccard_distance_matrix.qza")

AdNewt_BC16Smx<-read_qza("./16S/WithPlasmid/CoreMet_CopCorTrueAbundEst_r5000_Phylo/bray_curtis_distance_matrix.qza")

AdNewt_unWUF16Smx<-read_qza("./16S/WithPlasmid/CoreMet_CopCorTrueAbundEst_r5000_Phylo/unweighted_unifrac_distance_matrix.qza")

AdNewt_WUF16Smx<-read_qza("./16S/WithPlasmid/CoreMet_CopCorTrueAbundEst_r5000_Phylo/weighted_unifrac_distance_matrix.qza")

Merged Metadata, Alpha diversity and Antifungal predictions

AdNewt16S_Meta_Alpha_AntiF = Metadata16S %>%
  right_join(.,AdNewt_sOTU16S, by = "SampleID") %>%
  left_join(.,AdNewt_Evenness16S, by = "SampleID") %>%
  left_join(.,AdNewt_PD16S, by = "SampleID") %>%
  left_join(.,AdNewt_Shannon16S, by = "SampleID") %>%
  left_join(.,AntiFungal_Bsal, by = "SampleID") %>%
  left_join(.,AntiFungal_Full, by = "SampleID") %>%
  left_join(.,ReadCount16S, by = "SampleID") %>%
  mutate(Prop_AntiFungal_Count = TotalAntiFungal/CorrReadTotal)%>%
  mutate(Prop_AntiFungal_Bsal_Count = TotalAntiFungal_Bsal/CorrReadTotal) %>%
  mutate(ExpShannon = exp(shannon_entropy))

Alpha Diversity

Pre/Post Exposure

sOTU richness

#**Data wrangling to calculate change in alpha diversity through time **

AdNewt16S_Meta_Alpha_AntiF_PrePost = AdNewt16S_Meta_Alpha_AntiF %>%
  unite(DosInd, Dose, IndividualID, sep = "_", remove = FALSE) %>%
  select(SampleID,Temperature,TimeWeekCat,Dose,DosInd,observed_features,faith_pd,shannon_entropy,pielou_evenness, ExpShannon,AntiFungal_Richness,AntiFungal_Richness_Bsal,TotalAntiFungal,TotalAntiFungal_Bsal,Prop_AntiFungal_Count,Prop_AntiFungal_Bsal_Count) %>%
  filter(TimeWeekCat == "A" |TimeWeekCat == "B")

AdNewt16S_Meta_Alpha_AntiF_PrePost_6CRich = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "6") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,observed_features) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = observed_features) %>%
  mutate(ChangeAlpha = B-A)

AdNewt16S_Meta_Alpha_AntiF_PrePost_14CRich = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "14") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,observed_features) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = observed_features) %>%
  mutate(ChangeAlpha = B-A)

AdNewt16S_Meta_Alpha_AntiF_PrePost_22CRich = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "22") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,observed_features) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = observed_features) %>%
  mutate(ChangeAlpha = B-A)

PrePost_Rich_Differential<- rbind(AdNewt16S_Meta_Alpha_AntiF_PrePost_6CRich,AdNewt16S_Meta_Alpha_AntiF_PrePost_14CRich,AdNewt16S_Meta_Alpha_AntiF_PrePost_22CRich)

Dose <- c("Control","5x10.3","5x10.4","5x10.5","5x10.6")

DoseNum <- c("0","3","4","5","6")

DoseNumData <- data.frame(Dose,DoseNum)

PrePost_Rich_Differential$Temperature.F = factor(PrePost_Rich_Differential$Temperature, levels=c('6','14','22'))

PrePost_Rich_Differential %>%
  left_join(DoseNumData, by = "Dose") %>%
  ggplot(aes(x = as.numeric(DoseNum), y = ChangeAlpha)) +
  geom_point(aes(color = Temperature.F))+
  facet_grid(~Temperature.F) +
  theme(axis.text.x = element_text(angle = 90)) +
  xlab("Treatment Dose (10^x)")+
  ylab("Change in Richness (t1-t0) from pre to post infection")+
  geom_smooth(method=lm, level = 0.95, color = "black") +
  stat_cor(method = "kendall",label.x = 0.25, label.y = 80, label.sep = '\n')+
  scale_color_manual(values = c("#087CAF","#E7B800","#DB6E0B"))

Phylogenetic Diversity

AdNewt16S_Meta_Alpha_AntiF_PrePost_6Cfaith_pd = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "6") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,faith_pd) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = faith_pd) %>%
  mutate(ChangeAlpha = B-A)

AdNewt16S_Meta_Alpha_AntiF_PrePost_14Cfaith_pd = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "14") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,faith_pd) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = faith_pd) %>%
  mutate(ChangeAlpha = B-A)

AdNewt16S_Meta_Alpha_AntiF_PrePost_22Cfaith_pd = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "22") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,faith_pd) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = faith_pd) %>%
  mutate(ChangeAlpha = B-A)

PrePost_faith_pd_Differential<- rbind(AdNewt16S_Meta_Alpha_AntiF_PrePost_6Cfaith_pd,AdNewt16S_Meta_Alpha_AntiF_PrePost_14Cfaith_pd,AdNewt16S_Meta_Alpha_AntiF_PrePost_22Cfaith_pd)

PrePost_faith_pd_Differential$Temperature.F = factor(PrePost_faith_pd_Differential$Temperature, levels=c('6','14','22'))

PrePost_faith_pd_Differential %>%
  left_join(DoseNumData, by = "Dose") %>%
  ggplot(aes(x = as.numeric(DoseNum), y = ChangeAlpha)) +
  geom_point(aes(color = Temperature.F))+
  facet_grid(~Temperature.F) +
  theme(axis.text.x = element_text(angle = 90)) +
  xlab("Treatment Dose (10^x)")+
  ylab("Change in Faith's PD (t1-t0) from pre to post infection")+
  geom_smooth(method=lm, level = 0.95, color = "black") +
  stat_cor(method = "kendall",label.x = 0.25, label.sep = '\n')+
  #scale_colour_viridis_c(option = "plasma")+
  scale_color_manual(values = c("#087CAF","#E7B800","#DB6E0B"))

Shannon

AdNewt16S_Meta_Alpha_AntiF_PrePost_6Cshannon = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "6") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,shannon_entropy) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = shannon_entropy) %>%
  mutate(ChangeAlpha = B-A)

AdNewt16S_Meta_Alpha_AntiF_PrePost_14Cshannon = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "14") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,shannon_entropy) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = shannon_entropy) %>%
  mutate(ChangeAlpha = B-A)

AdNewt16S_Meta_Alpha_AntiF_PrePost_22Cshannon = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "22") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,shannon_entropy) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = shannon_entropy) %>%
  mutate(ChangeAlpha = B-A)

PrePost_shannon_Differential<- rbind(AdNewt16S_Meta_Alpha_AntiF_PrePost_6Cshannon,AdNewt16S_Meta_Alpha_AntiF_PrePost_14Cshannon,AdNewt16S_Meta_Alpha_AntiF_PrePost_22Cshannon)

PrePost_shannon_Differential$Temperature.F = factor(PrePost_shannon_Differential$Temperature, levels=c('6','14','22'))

PrePost_shannon_Differential %>%
  left_join(DoseNumData, by = "Dose") %>%
  ggplot(aes(x = as.numeric(DoseNum), y = ChangeAlpha)) +
geom_point(aes(color = Temperature.F))+
  facet_grid(~Temperature.F) +
  theme(axis.text.x = element_text(angle = 90)) +
  xlab("Treatment Dose (10^x)")+
  ylab("Change in Shannon (t1-t0) from pre to post infection")+
  geom_smooth(method=lm, level = 0.95, color = "black") +
  stat_cor(method = "kendall",label.x = 0.25, label.y = -2.5, label.sep = '\n')+
  #scale_colour_viridis_c(option = "plasma")+
  scale_color_manual(values = c("#087CAF","#E7B800","#DB6E0B"))

Effective species (exp-shannon)

AdNewt16S_Meta_Alpha_AntiF_PrePost_6CEffSp = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "6") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,ExpShannon) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = ExpShannon) %>%
  mutate(ChangeAlpha = B-A)

AdNewt16S_Meta_Alpha_AntiF_PrePost_14CEffSp = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "14") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,ExpShannon) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = ExpShannon) %>%
  mutate(ChangeAlpha = B-A)

AdNewt16S_Meta_Alpha_AntiF_PrePost_22CEffSp = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "22") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,ExpShannon) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = ExpShannon) %>%
  mutate(ChangeAlpha = B-A)

PrePost_EffSp_Differential<- rbind(AdNewt16S_Meta_Alpha_AntiF_PrePost_6CEffSp,AdNewt16S_Meta_Alpha_AntiF_PrePost_14CEffSp,AdNewt16S_Meta_Alpha_AntiF_PrePost_22CEffSp)

PrePost_EffSp_Differential$Temperature.F = factor(PrePost_EffSp_Differential$Temperature, levels=c('6','14','22'))

PrePost_EffSp_Differential %>%
  left_join(DoseNumData, by = "Dose") %>%
  ggplot(aes(x = as.numeric(DoseNum), y = ChangeAlpha)) +
  geom_point(aes(color = Temperature.F))+
  facet_grid(~Temperature.F) +
  theme(axis.text.x = element_text(angle = 90)) +
  xlab("Treatment Dose (10^x)")+
  ylab("Change in Effective Species (t1-t0) from pre to post infection")+
  geom_smooth(method=lm, level = 0.95, color = "black") +
  stat_cor(method = "kendall",label.x = 0.25, label.sep = '\n')+
  #scale_colour_viridis_c(option = "plasma")+
   scale_color_manual(values = c("#087CAF","#E7B800","#DB6E0B"))

Evenness

AdNewt16S_Meta_Alpha_AntiF_PrePost_6Cevenness = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "6") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,pielou_evenness) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = pielou_evenness) %>%
  mutate(ChangeAlpha = B-A)

AdNewt16S_Meta_Alpha_AntiF_PrePost_14Cevenness = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "14") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,pielou_evenness) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = pielou_evenness) %>%
  mutate(ChangeAlpha = B-A)

AdNewt16S_Meta_Alpha_AntiF_PrePost_22Cevenness = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "22") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,pielou_evenness) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = pielou_evenness) %>%
  mutate(ChangeAlpha = B-A)

PrePost_evenness_Differential<- rbind(AdNewt16S_Meta_Alpha_AntiF_PrePost_6Cevenness,AdNewt16S_Meta_Alpha_AntiF_PrePost_14Cevenness,AdNewt16S_Meta_Alpha_AntiF_PrePost_22Cevenness)

PrePost_evenness_Differential$Temperature.F = factor(PrePost_evenness_Differential$Temperature, levels=c('6','14','22'))

PrePost_evenness_Differential %>%
  left_join(DoseNumData, by = "Dose") %>%
  ggplot(aes(x = as.numeric(DoseNum), y = ChangeAlpha)) +
  geom_point(aes(color = Temperature.F))+
  facet_grid(~Temperature.F) +
  theme(axis.text.x = element_text(angle = 90)) +
  xlab("Treatment Dose (10^x)")+
  ylab("Change in Evenness (t1-t0) from pre to post infection")+
  geom_smooth(method=lm, level = 0.95, color = "black") +
  stat_cor(method = "kendall",label.x = 0.25, label.sep = '\n')+
  #scale_colour_viridis_c(option = "plasma")+
   scale_color_manual(values = c("#087CAF","#E7B800","#DB6E0B"))

Conclusions:

Overall: temperature induced differential responses of skin bacterial alpha diversity after exposure to Bsal

• at 14 C an increasing number of sOTUs and PD are lost as exposure dose increases

• at 6 C this pattern is reversed - increasing numbers of sOTUs and PD are gained with increasing exposure dose.

• at 22 C the change in richness is not correlated with dose for richness but following 14 C for PD

Full Time Series

sOTU Richness

size-scaled by Log Bsal load

AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Dose == "5x10.3" | Dose == "Control") %>%
  filter(!TimeWeek == "9")  %>%
  filter(!TimeWeek == "10")  %>%
  filter(!TimeWeek == "11") %>%
  #filter(!TimeWeek == "7") %>%
  #filter(!TimeWeek == "8") %>%
  #filter(!TimeWeek == "6") %>%
  #filter(!TimeWeek == "5") %>%
  ggplot(aes(TimeWeek,observed_features, color = Dose)) +
  geom_point(aes(size = LogCopies_uL))+
  #stat_summary(geom="errorbar", fun.data=mean_se, width=0) +
  #stat_summary(geom="point", fun.data=mean_se) +
  stat_summary(geom="line", fun.data=mean_se) +
  xlab("Week") +
  ylab(" Richness") +
  theme_q2r() + 
  scale_color_manual(values = c("#D1780F","#3F3F3E"))+
  facet_wrap(~Temp)

AdNewt16S_Meta_Alpha_AntiF_t1_4 = AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Dose == "5x10.3" | Dose == "Control") %>%
  filter(!TimeWeek == "9")  %>%
  filter(!TimeWeek == "10")  %>%
  filter(!TimeWeek == "11") %>%
  filter(!TimeWeek == "7") %>%
  filter(!TimeWeek == "8") %>%
  filter(!TimeWeek == "6") %>%
  filter(!TimeWeek == "5")

hist(AdNewt16S_Meta_Alpha_AntiF_t1_4$observed_features)

modsOTU_16St1_4 = lmer(observed_features ~ Dose*TimeWeek +Temperature + (1 | IndividualID), data = AdNewt16S_Meta_Alpha_AntiF_t1_4)

#sjPlot::plot_model(modsOTU_16St1_4, type = "diag")

car::Anova(modsOTU_16St1_4)

Conclusion: : Main Effect of Exposure, Temperature and Time

Contemplation: Model each temperature separately?

sOTU Richness

Only exposed treatment post exposure; size-scaled by Log Bsal load

AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Dose == "5x10.3"|Dose == "5x10.4"|Dose == "5x10.5"|Dose == "5x10.6") %>%
  filter(!TimeWeek == "0")  %>%
  filter(!TimeWeek == "9")  %>%
  filter(!TimeWeek == "10")  %>%
  filter(!TimeWeek == "11") %>%
  ggplot(aes(LogCopies_uL,observed_features)) +
  geom_point(aes(size = LogCopies_uL,color = Dose))+
  geom_smooth(method=lm, level = 0.95, color = "black") +
  stat_cor(method = "kendall",label.x = 1, label.y = 130,label.sep='\n')+
 #stat_summary(geom="errorbar", fun.data=mean_se, width=0) +
  #stat_summary(geom="point", fun.data=mean_se) +
  #stat_summary(geom="line", fun.data=mean_se) +
  xlab("Bsal") +
  ylab("Richness") +
  theme_q2r() + 
  facet_wrap(~Temp)

Phylogenetic Diversity

size-scaled by Log Bsal load

AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Dose == "5x10.3" | Dose == "Control") %>%
  filter(!TimeWeek == "9")  %>%
  filter(!TimeWeek == "10")  %>%
  filter(!TimeWeek == "11") %>%
  ggplot(aes(TimeWeek,faith_pd, color = Dose)) +
  geom_point(aes(size = LogCopies_uL))+
  #stat_summary(geom="errorbar", fun.data=mean_se, width=0) +
  #stat_summary(geom="point", fun.data=mean_se) +
  stat_summary(geom="line", fun.data=mean_se) +
  xlab("Week") +
  ylab(" Phylo Diversity") +
  theme_q2r() + 
  scale_color_manual(values = c("#D1780F","#3F3F3E"))+
  facet_wrap(~Temp)

Phylogenetic Diversity

Only exposed treatment; size-scaled by Log Bsal load

AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Dose == "5x10.3"|Dose == "5x10.4"|Dose == "5x10.5"|Dose == "5x10.6") %>%
  filter(!TimeWeek == "0")  %>%
  filter(!TimeWeek == "9")  %>%
  filter(!TimeWeek == "10")  %>%
  filter(!TimeWeek == "11") %>%
  ggplot(aes(LogCopies_uL,faith_pd)) +
  geom_point(aes(size = LogCopies_uL,color = Dose))+
  geom_smooth(method=lm, level = 0.95, color = "black") +
  stat_cor(method = "kendall",label.x = 1, label.y = 12,label.sep='\n')+
 #stat_summary(geom="errorbar", fun.data=mean_se, width=0) +
  #stat_summary(geom="point", fun.data=mean_se) +
  #stat_summary(geom="line", fun.data=mean_se) +
  xlab("Bsal") +
  ylab("Phylo Diversity") +
  theme_q2r() + 
  facet_wrap(~Temp)

Shannon Diversity

size-scaled by Log Bsal load

AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Dose == "5x10.3" | Dose == "Control") %>%
  filter(!TimeWeek == "9")  %>%
  filter(!TimeWeek == "10")  %>%
  filter(!TimeWeek == "11") %>%
  ggplot(aes(TimeWeek,shannon_entropy, color = Dose)) +
  geom_point(aes(size = LogCopies_uL))+
 #stat_summary(geom="errorbar", fun.data=mean_se, width=0) +
  #stat_summary(geom="point", fun.data=mean_se) +
  stat_summary(geom="line", fun.data=mean_se) +
  xlab("Week") +
  ylab("Shannon") +
  theme_q2r() + 
  scale_color_manual(values = c("#D1780F","#3F3F3E"))+
  facet_wrap(~Temp)

Shannon Diversity

Only exposed treatment; size-scaled by Log Bsal load

AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Dose == "5x10.3"|Dose == "5x10.4"|Dose == "5x10.5"|Dose == "5x10.6") %>%
  filter(!TimeWeek == "0")  %>%
  filter(!TimeWeek == "9")  %>%
  filter(!TimeWeek == "10")  %>%
  filter(!TimeWeek == "11") %>%
  ggplot(aes(LogCopies_uL,shannon_entropy)) +
  geom_point(aes(size = LogCopies_uL,color = Dose,))+
  geom_smooth(method=lm, level = 0.95, color = "black") +
  stat_cor(method = "kendall",label.x = 1, label.y = 5,label.sep='\n')+
 #stat_summary(geom="errorbar", fun.data=mean_se, width=0) +
  #stat_summary(geom="point", fun.data=mean_se) +
  #stat_summary(geom="line", fun.data=mean_se) +
  xlab("Bsal") +
  ylab("Shannon") +
  theme_q2r() + 
  facet_wrap(~Temp)

Evenness

size-scaled by Log Bsal load

AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Dose == "5x10.3" | Dose == "Control") %>%
  filter(!TimeWeek == "9")  %>%
  filter(!TimeWeek == "10")  %>%
  filter(!TimeWeek == "11") %>%
  ggplot(aes(TimeWeek,pielou_evenness, color = Dose)) +
  geom_point(aes(size = LogCopies_uL))+
  #stat_summary(geom="errorbar", fun.data=mean_se, width=0) +
  #stat_summary(geom="point", fun.data=mean_se) +
  stat_summary(geom="line", fun.data=mean_se) +
  xlab("Days") +
  ylab(" Evenness") +
  theme_q2r() + 
  scale_color_manual(values = c("#D1780F","#3F3F3E"))+
  facet_wrap(~Temp)

Evenness

Only exposed treatment; size-scaled by Log Bsal load

AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Dose == "5x10.3"|Dose == "5x10.4"|Dose == "5x10.5"|Dose == "5x10.6") %>%
  filter(!TimeWeek == "0")  %>%
  filter(!TimeWeek == "9")  %>%
  filter(!TimeWeek == "10")  %>%
  filter(!TimeWeek == "11") %>%
  ggplot(aes(LogCopies_uL,pielou_evenness)) +
  geom_point(aes(size = LogCopies_uL,color = Dose))+
  geom_smooth(method=lm, level = 0.95, color = "black") +
  stat_cor(method = "kendall",label.x = 1, label.y = .25,label.sep='\n')+
 #stat_summary(geom="errorbar", fun.data=mean_se, width=0) +
  #stat_summary(geom="point", fun.data=mean_se) +
  #stat_summary(geom="line", fun.data=mean_se) +
  xlab("Bsal") +
  ylab("Shannon") +
  theme_q2r() + 
  facet_wrap(~Temp)

Predicted Function

Pre/Post Exposure

Antifungal Richness

AdNewt16S_Meta_Alpha_AntiF_PrePost_6C_AFRich = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "6") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,AntiFungal_Richness) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = AntiFungal_Richness) %>%
  mutate(ChangeAlpha = B-A)

AdNewt16S_Meta_Alpha_AntiF_PrePost_14C_AFRich = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "14") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,AntiFungal_Richness) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = AntiFungal_Richness) %>%
  mutate(ChangeAlpha = B-A)

AdNewt16S_Meta_Alpha_AntiF_PrePost_22C_AFRich = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "22") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,AntiFungal_Richness) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = AntiFungal_Richness) %>%
  mutate(ChangeAlpha = B-A)

PrePost_AFRich_Differential<- rbind(AdNewt16S_Meta_Alpha_AntiF_PrePost_6C_AFRich,AdNewt16S_Meta_Alpha_AntiF_PrePost_14C_AFRich,AdNewt16S_Meta_Alpha_AntiF_PrePost_22C_AFRich)

PrePost_AFRich_Differential$Temperature.F = factor(PrePost_AFRich_Differential$Temperature, levels=c('6','14','22'))

PrePost_AFRich_Differential %>%
  left_join(DoseNumData, by = "Dose") %>%
  ggplot(aes(x = as.numeric(DoseNum), y = ChangeAlpha)) +
  geom_point(aes(color = Temperature.F))+
  facet_grid(~Temperature.F) +
  theme(axis.text.x = element_text(angle = 90)) +
  xlab("Treatment Dose (10^x)")+
  ylab("Change in Antifungal Richness (t1-t0) from pre to post infection")+
  geom_smooth(method=lm, level = 0.95, color = "black") +
  stat_cor(method = "kendall",label.x = 0.25, label.sep = '\n')+
  #scale_colour_viridis_c(option = "plasma")+
  scale_color_manual(values = c("#087CAF","#E7B800","#DB6E0B"))

Antifungal Function from corrected unrarified data

AdNewt16S_Meta_Alpha_AntiF_PrePost_6CTotalAntiFungal = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "6") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,TotalAntiFungal) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = TotalAntiFungal) %>%
  mutate(ChangeAlpha = B-A)

AdNewt16S_Meta_Alpha_AntiF_PrePost_14CTotalAntiFungal = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "14") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,TotalAntiFungal) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = TotalAntiFungal) %>%
  mutate(ChangeAlpha = B-A)

AdNewt16S_Meta_Alpha_AntiF_PrePost_22CTotalAntiFungal = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "22") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,TotalAntiFungal) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = TotalAntiFungal) %>%
  mutate(ChangeAlpha = B-A)

PrePost_TotalAntiFungal_Differential<- rbind(AdNewt16S_Meta_Alpha_AntiF_PrePost_6CTotalAntiFungal,AdNewt16S_Meta_Alpha_AntiF_PrePost_14CTotalAntiFungal,AdNewt16S_Meta_Alpha_AntiF_PrePost_22CTotalAntiFungal)

PrePost_TotalAntiFungal_Differential$Temperature.F = factor(PrePost_TotalAntiFungal_Differential$Temperature, levels=c('6','14','22'))

PrePost_TotalAntiFungal_Differential %>%
  left_join(DoseNumData, by = "Dose") %>%
  ggplot(aes(x = as.numeric(DoseNum), y = log10(ChangeAlpha))) +
  geom_point(aes(color = Temperature.F))+
  facet_grid(~Temperature.F) +
  theme(axis.text.x = element_text(angle = 90)) +
  xlab("Treatment Dose (10^x)")+
  ylab("Change in Antifungal Function (t1-t0) from pre to post infection")+
  geom_smooth(method=lm, level = 0.95, color = "black") +
  stat_cor(method = "kendall",label.x = 0.25, label.sep = '\n')+
  #scale_colour_viridis_c(option = "plasma")+
  scale_color_manual(values = c("#087CAF","#E7B800","#DB6E0B"))

Antifungal Function

proportionalized through total corrected “abundance” estimate

AdNewt16S_Meta_Alpha_AntiF_PrePost_6CProp_AntiFungal_Count = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "6") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,Prop_AntiFungal_Count) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = Prop_AntiFungal_Count) %>%
  mutate(ChangeAlpha = B-A)

AdNewt16S_Meta_Alpha_AntiF_PrePost_14CProp_AntiFungal_Count = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "14") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,Prop_AntiFungal_Count) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = Prop_AntiFungal_Count) %>%
  mutate(ChangeAlpha = B-A)

AdNewt16S_Meta_Alpha_AntiF_PrePost_22CProp_AntiFungal_Count = AdNewt16S_Meta_Alpha_AntiF_PrePost %>%
 filter(Temperature == "22") %>%
  select(Temperature,Dose,TimeWeekCat,DosInd,Prop_AntiFungal_Count) %>%
  pivot_wider(names_from = TimeWeekCat,values_from = Prop_AntiFungal_Count) %>%
  mutate(ChangeAlpha = B-A)

PrePost_TotalAntiFungalCt_Differential<- rbind(AdNewt16S_Meta_Alpha_AntiF_PrePost_6CProp_AntiFungal_Count,AdNewt16S_Meta_Alpha_AntiF_PrePost_14CProp_AntiFungal_Count,AdNewt16S_Meta_Alpha_AntiF_PrePost_22CProp_AntiFungal_Count)

PrePost_TotalAntiFungalCt_Differential$Temperature.F = factor(PrePost_TotalAntiFungalCt_Differential$Temperature, levels=c('6','14','22'))

PrePost_TotalAntiFungalCt_Differential %>%
  left_join(DoseNumData, by = "Dose") %>%
  ggplot(aes(x = as.numeric(DoseNum), y = ChangeAlpha)) +
  geom_point(aes(color = Temperature.F))+
  facet_grid(~Temperature.F) +
  theme(axis.text.x = element_text(angle = 90)) +
  xlab("Treatment Dose (10^x)")+
  ylab("Change in Proportional Function (t1-t0) from pre to post infection")+
  geom_smooth(method=lm, level = 0.95, color = "black") +
  stat_cor(method = "pearson",label.x = 0.25,label.y = -0.25, label.sep = '\n')+
  #scale_colour_viridis_c(option = "plasma")+
   scale_color_manual(values = c("#087CAF","#E7B800","#DB6E0B"))

Antifungal function of between pre/post exposure: 6 C

AdNewt16S_Meta_Alpha_AntiF %>%
  #filter(Dose=="Control") %>%
  #filter(TimeWeek == "0") %>%
  filter(TimeWeek == "0"|TimeWeek=="1") %>%
  filter(Temp=="6")%>%
  ggplot(aes(TimeWeek,TotalAntiFungal, color = Dose,)) +
  geom_point(aes(size = LogCopies_uL,alpha = 0.5))+
  ylab("Antifungal Function") +
  theme_q2r()+
  geom_smooth(method=lm, level = 0.95) +
  stat_cor(method = "kendall",label.x = 0, label.y = 0.8, label.sep='\n')+
  facet_grid(~factor(Dose, levels=c('Control','5x10.3','5x10.4','5x10.5', '5x10.6')))+
  ggtitle("6 C Newts")

Antifungal function of between pre/post exposure: 14 C

AdNewt16S_Meta_Alpha_AntiF %>%
  #filter(Dose=="Control") %>%
  #filter(TimeWeek == "0") %>%
  filter(TimeWeek == "0"|TimeWeek=="1") %>%
  filter(Temp=="14")%>%
  ggplot(aes(TimeWeek,TotalAntiFungal, color = Dose)) +
  geom_point(aes(size = LogCopies_uL,alpha = 0.5))+
  ylab("Antifungal Function") +
  theme_q2r()+
  geom_smooth(method=lm, level = 0.95) +
  stat_cor(method = "kendall",label.x = 0, label.y = 0.8, label.sep='\n')+
  facet_grid(~factor(Dose, levels=c('Control','5x10.3','5x10.4','5x10.5', '5x10.6')))+
  ggtitle("14 C Newts")

Antifungal function of between pre/post exposure: 22 C

AdNewt16S_Meta_Alpha_AntiF %>%
  #filter(Dose=="Control") %>%
  #filter(TimeWeek == "0") %>%
  filter(TimeWeek == "0"|TimeWeek=="1") %>%
  filter(Temp=="22")%>%
  ggplot(aes(TimeWeek,TotalAntiFungal, color = Dose)) +
  geom_point()+
  ylab("Antifungal Function") +
  theme_q2r()+
  geom_smooth(method=lm, level = 0.95) +
  stat_cor(method = "kendall",label.x = 0, label.y = 0.8, label.sep='\n')+
  facet_grid(~factor(Dose, levels=c('Control','5x10.3','5x10.4','5x10.5', '5x10.6')))+
  ggtitle("22 C Newts")

Full Time Series

Antifungal Function ~ total corrected abundance

AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Dose == "5x10.3" | Dose == "Control") %>%
  filter(TimeWeek<8)  %>%
  filter(!TempTime=="14C.6") %>%
  filter(!TempTime=="14C.7") %>%
  ggplot(aes(TimeWeek,log10(TotalAntiFungal), color = Dose)) +
  stat_summary(geom="errorbar", fun.data=mean_se, width=0) +
  stat_summary(geom="point", fun.data=mean_se) +
  stat_summary(geom="line", fun.data=mean_se) +
  xlab("Days") +
  ylab(" log Antifungal Function") +
  theme_q2r() + 
  scale_color_colorblind()+
  facet_wrap(~Temp)

Antifungal Function ~ total corrected abundance

Only exposed treatment; size-scaled by Log Bsal load

AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Dose == "5x10.3"|Dose == "5x10.4"|Dose == "5x10.5"|Dose == "5x10.6") %>%
  filter(!TimeWeek == "0")  %>%
  filter(!TimeWeek == "9")  %>%
  filter(!TimeWeek == "10")  %>%
  filter(!TimeWeek == "11") %>%
  ggplot(aes(LogCopies_uL,log10(TotalAntiFungal))) +
  geom_point(aes(size = LogCopies_uL,color = Dose))+
  geom_smooth(method=lm, level = 0.95, color = "black") +
  stat_cor(method = "kendall",label.x = 1, label.y = 3,label.sep='\n')+
 #stat_summary(geom="errorbar", fun.data=mean_se, width=0) +
  #stat_summary(geom="point", fun.data=mean_se) +
  #stat_summary(geom="line", fun.data=mean_se) +
  xlab("Bsal") +
  ylab("Antifungal Function") +
  theme_q2r() + 
  facet_wrap(~Temp)

Beta Diversity

Location

AdNewt16S_Meta_Alpha_AntiF_T0 <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(TimeWeekCat== "A") %>%
  drop_na(observed_features)
  
metaList_T0<-as.vector(AdNewt16S_Meta_Alpha_AntiF_T0$SampleID)

AdNewt_unWUF16Smx2 = AdNewt_unWUF16Smx$data

AdNewt_unWUF16Smx_T0 <- usedist::dist_subset(AdNewt_unWUF16Smx2,metaList_T0)

AdNewt_unWUF16Smx_T0b <-as.matrix(AdNewt_unWUF16Smx_T0)

PCoA_AdNewt_unWUF16Smx_T0b <- ape::pcoa(AdNewt_unWUF16Smx_T0b)

PCoA_AdNewt_unWUF16Smx_T0b_Axes = as.data.frame(PCoA_AdNewt_unWUF16Smx_T0b$vectors)
PCoA_AdNewt_unWUF16Smx_T0b_Axesb =tibble::rownames_to_column(PCoA_AdNewt_unWUF16Smx_T0b_Axes, "SampleID")

PCoA_AdNewt_unWUF16Smx_T0b_Axesb %>%
  select(SampleID, Axis.1, Axis.2) %>%
  left_join(AdNewt16S_Meta_Alpha_AntiF, by = "SampleID") %>%
  ggplot(aes(x=Axis.1, y=Axis.2)) +
  geom_point(aes(color=CollectionLocation, shape = Temperature, alpha=0.5,size =2))+
  theme_q2r() +
  scale_color_calc(name="Location") +
  #stat_ellipse(type = "t",linetype = 2)+
  ggtitle("All temps - unweighted")

PCoA_AdNewt_unWUF16Smx_T0b_Axesb %>%
  select(SampleID, Axis.1, Axis.2) %>%
  left_join(AdNewt16S_Meta_Alpha_AntiF, by = "SampleID") %>%
  ggplot(aes(x=Axis.1, y=Axis.2)) +
  geom_point(aes(color=CollectionLocation, alpha=0.5,size =2))+
  theme_q2r() +
  scale_color_calc(name="Location") +
  #stat_ellipse(type = "t",linetype = 2)+
  ggtitle("All temps - unweighted")+
  facet_wrap(~Temperature)

Pre/Post

Temperature 6 C

AdNewt16S_Meta_Alpha_AntiF_6_T0_T1 <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 6) %>%
  filter(TimeWeekCat== "A"|TimeWeekCat=="B") %>%
  drop_na(observed_features)
  
metaList_6t0t1<-as.vector(AdNewt16S_Meta_Alpha_AntiF_6_T0_T1$SampleID)

AdNewt_unWUF16Smx2 = AdNewt_unWUF16Smx$data

AdNewt_unWUF16Smx_6t0t1 <- usedist::dist_subset(AdNewt_unWUF16Smx2,metaList_6t0t1)

AdNewt_unWUF16Smx_6t0t1b <-as.matrix(AdNewt_unWUF16Smx_6t0t1)

AdNewt6t0t1_16S_distList<-colnames(AdNewt_unWUF16Smx_6t0t1b)

AdNewt16S_Meta_Alpha_AntiF_6_T0_T1b<- AdNewt16S_Meta_Alpha_AntiF_6_T0_T1 %>%
  filter(SampleID %in% AdNewt6t0t1_16S_distList)

PCoA_AdNewt6t0t1_16S <- ape::pcoa(AdNewt_unWUF16Smx_6t0t1)

PCoA_AdNewt6t0t1_16S_Axes = as.data.frame(PCoA_AdNewt6t0t1_16S$vectors)
PCoA_AdNewt6t0t1_16S_Axesb =tibble::rownames_to_column(PCoA_AdNewt6t0t1_16S_Axes, "SampleID")

PCoA_AdNewt6t0t1_16S_Axesb %>%
  select(SampleID, Axis.1, Axis.2) %>%
  left_join(AdNewt16S_Meta_Alpha_AntiF, by = "SampleID") %>%
  ggplot(aes(x=Axis.1, y=Axis.2)) +
  geom_point(aes(color=Dose, alpha=0.5,shape =TimeWeekCat,size =2))+
  theme_q2r() +
  scale_color_colorblind(name="TimeWeek") +
  #stat_ellipse(type = "t",linetype = 2)+
  ggtitle("6C - unweighted")

AdNewt16S_Meta_Alpha_AntiF$Dose = factor(AdNewt16S_Meta_Alpha_AntiF$Dose, levels=c('Control','5x10.3','5x10.4','5x10.5','5x10.6'))

PCoA_AdNewt6t0t1_16S_Axesb %>%
  select(SampleID, Axis.1, Axis.2) %>%
  left_join(AdNewt16S_Meta_Alpha_AntiF, by = "SampleID") %>%
  ggplot(aes(x=Axis.1, y=Axis.2,shape =TimeWeekCat)) +
  geom_point(aes(color=Dose, alpha=0.5,shape =TimeWeekCat,size =2))+
  theme_q2r() +
  scale_color_colorblind(name="TimeWeek") +
  stat_ellipse(type = "t",linetype = 2)+
  ggtitle("6C - unweighted")+
  facet_wrap(~Dose)

Temperature 14 C

AdNewt16S_Meta_Alpha_AntiF_14_T0_T1 <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 14) %>%
  filter(TimeWeekCat== "A"|TimeWeekCat=="B") %>%
  drop_na(observed_features)
  
metaList_14t0t1<-as.vector(AdNewt16S_Meta_Alpha_AntiF_14_T0_T1$SampleID)

AdNewt_unWUF16Smx2 = AdNewt_unWUF16Smx$data

AdNewt_unWUF16Smx_14t0t1 <- usedist::dist_subset(AdNewt_unWUF16Smx2,metaList_14t0t1)

AdNewt_unWUF16Smx_14t0t1b <-as.matrix(AdNewt_unWUF16Smx_14t0t1)

AdNewt14t0t1_16S_distList<-colnames(AdNewt_unWUF16Smx_14t0t1b)

AdNewt16S_Meta_Alpha_AntiF_14_T0_T1b<- AdNewt16S_Meta_Alpha_AntiF_14_T0_T1 %>%
  filter(SampleID %in% AdNewt14t0t1_16S_distList)

PCoA_AdNewt14t0t1_16S <- ape::pcoa(AdNewt_unWUF16Smx_14t0t1)

PCoA_AdNewt14t0t1_16S_Axes = as.data.frame(PCoA_AdNewt14t0t1_16S$vectors)
PCoA_AdNewt14t0t1_16S_Axesb =tibble::rownames_to_column(PCoA_AdNewt14t0t1_16S_Axes, "SampleID")

PCoA_AdNewt14t0t1_16S_Axesb %>%
  select(SampleID, Axis.1, Axis.2) %>%
  left_join(AdNewt16S_Meta_Alpha_AntiF, by = "SampleID") %>%
  ggplot(aes(x=Axis.1, y=Axis.2)) +
  geom_point(aes(color=Dose, alpha=0.5,shape =TimeWeekCat,size =2))+
  theme_q2r() +
  scale_color_colorblind(name="TimeWeek") +
  #stat_ellipse(type = "t",linetype = 2)+
  ggtitle("14C - unweighted")

AdNewt16S_Meta_Alpha_AntiF$Dose = factor(AdNewt16S_Meta_Alpha_AntiF$Dose, levels=c('Control','5x10.3','5x10.4','5x10.5','5x10.6'))

PCoA_AdNewt14t0t1_16S_Axesb %>%
  select(SampleID, Axis.1, Axis.2) %>%
  left_join(AdNewt16S_Meta_Alpha_AntiF, by = "SampleID") %>%
  ggplot(aes(x=Axis.1, y=Axis.2,shape =TimeWeekCat)) +
  geom_point(aes(color=Dose, alpha=0.5,shape =TimeWeekCat,size =2))+
  theme_q2r() +
  scale_color_colorblind(name="TimeWeek") +
  stat_ellipse(type = "t",linetype = 2)+
  ggtitle("14C - unweighted")+
  facet_wrap(~Dose)

Temperature 22

AdNewt16S_Meta_Alpha_AntiF_22_T0_T1 <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 22) %>%
  filter(TimeWeekCat== "A"|TimeWeekCat=="B") %>%
  drop_na(observed_features)
  
metaList_22t0t1<-as.vector(AdNewt16S_Meta_Alpha_AntiF_22_T0_T1$SampleID)

AdNewt_unWUF16Smx2 = AdNewt_unWUF16Smx$data

AdNewt_unWUF16Smx_22t0t1 <- usedist::dist_subset(AdNewt_unWUF16Smx2,metaList_22t0t1)

AdNewt_unWUF16Smx_22t0t1b <-as.matrix(AdNewt_unWUF16Smx_22t0t1)

AdNewt22t0t1_16S_distList<-colnames(AdNewt_unWUF16Smx_22t0t1b)

#AdNewt16S_Meta_Alpha_AntiF_22_T0_T1b<- AdNewt16S_Meta_Alpha_AntiF_22_T0_T1 %>%
  #filter(SampleID %in% AdNewt22t0t1_16S_distList)

PCoA_AdNewt22t0t1_16S <- ape::pcoa(AdNewt_unWUF16Smx_22t0t1)

PCoA_AdNewt22t0t1_16S_Axes = as.data.frame(PCoA_AdNewt22t0t1_16S$vectors)
PCoA_AdNewt22t0t1_16S_Axesb =tibble::rownames_to_column(PCoA_AdNewt22t0t1_16S_Axes, "SampleID")

PCoA_AdNewt22t0t1_16S_Axesb %>%
  select(SampleID, Axis.1, Axis.2) %>%
  left_join(AdNewt16S_Meta_Alpha_AntiF, by = "SampleID") %>%
  ggplot(aes(x=Axis.1, y=Axis.2)) +
  geom_point(aes(color=Dose, alpha=0.5,shape =TimeWeekCat,size =2))+
  theme_q2r() +
  scale_color_colorblind(name="TimeWeek") +
  #stat_ellipse(type = "t",linetype = 2)+
  ggtitle("22C - unweighted")

PCoA_AdNewt22t0t1_16S_Axesb %>%
  select(SampleID, Axis.1, Axis.2) %>%
  left_join(AdNewt16S_Meta_Alpha_AntiF, by = "SampleID") %>%
  ggplot(aes(x=Axis.1, y=Axis.2,shape =TimeWeekCat)) +
  geom_point(aes(color=Dose, alpha=0.5,shape =TimeWeekCat,size =2))+
  theme_q2r() +
  scale_color_colorblind(name="TimeWeek") +
  stat_ellipse(type = "t",linetype = 2)+
  ggtitle("22C - unweighted")+
  facet_wrap(~Dose)

Meta16S_14_t1t2_Con <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 14) %>%
  filter(TimeWeekCat== "A"|TimeWeekCat=="B") %>%
  filter(Dose=="Control")%>%
  drop_na(observed_features)
  
metaList_14t1t2_Con<-as.vector(Meta16S_14_t1t2_Con$SampleID)

AdNewt_unWUF16Smx_14t0t1Con <- usedist::dist_subset(AdNewt_unWUF16Smx2,metaList_14t1t2_Con)
AdNewt_unWUF16Smx_14t0t1Conb <-as.matrix(AdNewt_unWUF16Smx_14t0t1Con)


AdNewt16S_Meta_Alpha_AntiF

PairwiseData_16S_T14Conab <- adegenet::pairDist(AdNewt_unWUF16Smx_14t0t1Conb,Meta16S_14_t1t2_Con$TimeWeekCat, within = FALSE)

PairwiseData_16S_T14Conab2 <- PairwiseData_16S_T14Conab$data %>%
  add_column(Temperature = "14")%>%
  add_column(Dose= "Control")

## 3 dose

Meta16S_14_t1t2_3 <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 14) %>%
  filter(TimeWeekCat== "A"|TimeWeekCat=="B") %>%
  filter(Dose=="5x10.3")%>%
  drop_na(observed_features)
  
metaList_14t1t2_3<-as.vector(Meta16S_14_t1t2_3$SampleID)

AdNewt_unWUF16Smx_14t0t13 <- usedist::dist_subset(AdNewt_unWUF16Smx2,metaList_14t1t2_3)
AdNewt_unWUF16Smx_14t0t13b <-as.matrix(AdNewt_unWUF16Smx_14t0t13)



PairwiseData_16S_T143ab <- adegenet::pairDist(AdNewt_unWUF16Smx_14t0t13b,Meta16S_14_t1t2_3$TimeWeekCat, within = FALSE)

PairwiseData_16S_T143ab2 <- PairwiseData_16S_T143ab$data %>%
  add_column(Temperature = "14")%>%
  add_column(Dose= "5x10.3")

## 4 
Meta16S_14_t1t2_4 <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 14) %>%
  filter(TimeWeekCat== "A"|TimeWeekCat=="B") %>%
  filter(Dose=="5x10.4")%>%
  drop_na(observed_features)
  
metaList_14t1t2_4<-as.vector(Meta16S_14_t1t2_4$SampleID)

AdNewt_unWUF16Smx_14t0t14 <- usedist::dist_subset(AdNewt_unWUF16Smx2,metaList_14t1t2_4)
AdNewt_unWUF16Smx_14t0t14b <-as.matrix(AdNewt_unWUF16Smx_14t0t14)

PairwiseData_16S_T144ab <- adegenet::pairDist(AdNewt_unWUF16Smx_14t0t14b,Meta16S_14_t1t2_4$TimeWeekCat, within = FALSE)

PairwiseData_16S_T144ab2 <- PairwiseData_16S_T144ab$data %>%
  add_column(Temperature = "14")%>%
  add_column(Dose= "5x10.4")

## 5 
Meta16S_14_t1t2_5 <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 14) %>%
  filter(TimeWeekCat== "A"|TimeWeekCat=="B") %>%
  filter(Dose=="5x10.5")%>%
  drop_na(observed_features)
  
metaList_14t1t2_5<-as.vector(Meta16S_14_t1t2_5$SampleID)

AdNewt_unWUF16Smx_14t0t15 <- usedist::dist_subset(AdNewt_unWUF16Smx2,metaList_14t1t2_5)
AdNewt_unWUF16Smx_14t0t15b <-as.matrix(AdNewt_unWUF16Smx_14t0t15)

PairwiseData_16S_T145ab <- adegenet::pairDist(AdNewt_unWUF16Smx_14t0t15b,Meta16S_14_t1t2_5$TimeWeekCat, within = FALSE)

PairwiseData_16S_T145ab2 <- PairwiseData_16S_T145ab$data %>%
  add_column(Temperature = "14")%>%
  add_column(Dose= "5x10.5")

## 6
Meta16S_14_t1t2_6 <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 14) %>%
  filter(TimeWeekCat== "A"|TimeWeekCat=="B") %>%
  filter(Dose=="5x10.6")%>%
  drop_na(observed_features)
  
metaList_14t1t2_6<-as.vector(Meta16S_14_t1t2_6$SampleID)

AdNewt_unWUF16Smx_14t0t16 <- usedist::dist_subset(AdNewt_unWUF16Smx2,metaList_14t1t2_6)
AdNewt_unWUF16Smx_14t0t16b <-as.matrix(AdNewt_unWUF16Smx_14t0t16)

PairwiseData_16S_T146ab <- adegenet::pairDist(AdNewt_unWUF16Smx_14t0t16b,Meta16S_14_t1t2_6$TimeWeekCat, within = FALSE)

PairwiseData_16S_T146ab2 <- PairwiseData_16S_T146ab$data %>%
  add_column(Temperature = "14")%>%
  add_column(Dose= "5x10.6")

Pairwise_14_Compile = rbind(PairwiseData_16S_T14Conab2,PairwiseData_16S_T143ab2,PairwiseData_16S_T144ab2,PairwiseData_16S_T145ab2,PairwiseData_16S_T146ab2)

Pairwise_14_Compile  %>%
  left_join(.,DoseNumData) %>%
  ggplot(aes(x = as.numeric(DoseNum),y = distance))+
  geom_point()+ 
  geom_smooth(method=lm, level = 0.95)+ stat_cor(method = "kendall", label.sep='\n') +
  ggtitle("Pairwise distance between t0 and t1 across doses")+
  xlab("Dose")

Time Series

Temperature 6 - through time

AdNewt16S_Meta_Alpha_AntiF_6_TS <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 6) %>%
  filter(Dose == "5x10.3" | Dose == "Control") %>%
  filter(TimeWeek<8) %>%
  drop_na(observed_features)
  
metaList_6_TS<-as.vector(AdNewt16S_Meta_Alpha_AntiF_6_TS$SampleID)

AdNewt_unWUF16Smx2 = AdNewt_unWUF16Smx$data

AdNewt_unWUF16Smx_6_TS <- usedist::dist_subset(AdNewt_unWUF16Smx2,metaList_6_TS)

AdNewt_unWUF16Smx_6_TSb <-as.matrix(AdNewt_unWUF16Smx_6_TS)

AdNewt6_TS_16S_distList<-colnames(AdNewt_unWUF16Smx_6_TSb)

AdNewt16S_Meta_Alpha_AntiF_6_TSb<- AdNewt16S_Meta_Alpha_AntiF_6_TS %>%
  filter(SampleID %in% AdNewt6_TS_16S_distList)

PCoA_AdNewt_unWUF16Smx_6_TSb <- ape::pcoa(AdNewt_unWUF16Smx_6_TSb)

PCoA_AdNewt_unWUF16Smx_6_TSb_Axes = as.data.frame(PCoA_AdNewt_unWUF16Smx_6_TSb$vectors)
PCoA_AdNewt_unWUF16Smx_6_TSb_Axesb =tibble::rownames_to_column(PCoA_AdNewt_unWUF16Smx_6_TSb_Axes, "SampleID")

PCoA_AdNewt_unWUF16Smx_6_TSb_Axesb %>%
  select(SampleID, Axis.1, Axis.2,Axis.3) %>%
  left_join(AdNewt16S_Meta_Alpha_AntiF, by = "SampleID") %>%
  ggplot(aes(x=Axis.1, y=Axis.2)) +
  geom_point(aes(color=TimeWeek, alpha=0.5,shape =Dose,size =LogCopies_uL))+
  theme_q2r() +
  scale_color_viridis_c(option = "viridis", name="TimeWeek") +
  #stat_ellipse(type = "t",linetype = 2)+
  ggtitle("6C - unweighted")

PCoA_AdNewt_unWUF16Smx_6_TSb_Axesb %>%
  select(SampleID, Axis.1, Axis.2,Axis.3) %>%
  left_join(AdNewt16S_Meta_Alpha_AntiF, by = "SampleID") %>%
  ggplot(aes(x=Axis.1, y=Axis.2,color=Dose,shape =Dose)) +
  geom_point(aes(alpha=0.8,shape =Dose,size =LogCopies_uL))+
  theme_q2r() +
  scale_color_colorblind() +
  #stat_ellipse(type = "t",linetype = 2)+
  ggtitle("6C - unweighted")+
  stat_ellipse(type = "t",linetype = 2) +
  facet_wrap(~TimeWeek)

Temperature 14 - through time

AdNewt16S_Meta_Alpha_AntiF_14_TS <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 14) %>%
  filter(Dose == "5x10.3" | Dose == "Control") %>%
  filter(TimeWeek<5) %>%
  drop_na(observed_features)
  
metaList_14_TS<-as.vector(AdNewt16S_Meta_Alpha_AntiF_14_TS$SampleID)

AdNewt_unWUF16Smx2 = AdNewt_unWUF16Smx$data

AdNewt_unWUF16Smx_14_TS <- usedist::dist_subset(AdNewt_unWUF16Smx2,metaList_14_TS)

AdNewt_unWUF16Smx_14_TSb <-as.matrix(AdNewt_unWUF16Smx_14_TS)

AdNewt14_TS_16S_distList<-colnames(AdNewt_unWUF16Smx_14_TSb)

AdNewt16S_Meta_Alpha_AntiF_14_TSb<- AdNewt16S_Meta_Alpha_AntiF_14_TS %>%
  filter(SampleID %in% AdNewt14_TS_16S_distList)

PCoA_AdNewt_unWUF16Smx_14_TSb <- ape::pcoa(AdNewt_unWUF16Smx_14_TSb)

PCoA_AdNewt_unWUF16Smx_14_TSb_Axes = as.data.frame(PCoA_AdNewt_unWUF16Smx_14_TSb$vectors)
PCoA_AdNewt_unWUF16Smx_14_TSb_Axesb =tibble::rownames_to_column(PCoA_AdNewt_unWUF16Smx_14_TSb_Axes, "SampleID")

PCoA_AdNewt_unWUF16Smx_14_TSb_Axesb %>%
  select(SampleID, Axis.1, Axis.2,Axis.3) %>%
  left_join(AdNewt16S_Meta_Alpha_AntiF, by = "SampleID") %>%
  ggplot(aes(x=Axis.1, y=Axis.2)) +
  geom_point(aes(color=TimeWeek, alpha=0.5,shape =Dose,size =LogCopies_uL))+
  theme_q2r() +
  scale_color_viridis_c(option = "viridis", name="TimeWeek") +
  #stat_ellipse(type = "t",linetype = 2)+
  ggtitle("14C - unweighted")

PCoA_AdNewt_unWUF16Smx_14_TSb_Axesb %>%
  select(SampleID, Axis.1, Axis.2,Axis.3) %>%
  left_join(AdNewt16S_Meta_Alpha_AntiF, by = "SampleID") %>%
  ggplot(aes(x=Axis.1, y=Axis.2,color=Dose,shape =Dose)) +
  geom_point(aes(alpha=0.8,shape =Dose,size =LogCopies_uL))+
  theme_q2r() +
  scale_color_colorblind() +
  #stat_ellipse(type = "t",linetype = 2)+
  ggtitle("14C - unweighted")+
  stat_ellipse(type = "t",linetype = 2) +
  facet_wrap(~TimeWeek)

Temperature 22 - through time

AdNewt16S_Meta_Alpha_AntiF_22_TS <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 22) %>%
  filter(Dose == "5x10.3" | Dose == "Control") %>%
  filter(TimeWeek<8) %>%
  drop_na(observed_features)
  
metaList_22_TS<-as.vector(AdNewt16S_Meta_Alpha_AntiF_22_TS$SampleID)

AdNewt_unWUF16Smx2 = AdNewt_unWUF16Smx$data

AdNewt_unWUF16Smx_22_TS <- usedist::dist_subset(AdNewt_unWUF16Smx2,metaList_22_TS)

AdNewt_unWUF16Smx_22_TSb <-as.matrix(AdNewt_unWUF16Smx_22_TS)

AdNewt22_TS_16S_distList<-colnames(AdNewt_unWUF16Smx_22_TSb)

AdNewt16S_Meta_Alpha_AntiF_22_TSb<- AdNewt16S_Meta_Alpha_AntiF_22_TS %>%
  filter(SampleID %in% AdNewt22_TS_16S_distList)

PCoA_AdNewt_unWUF16Smx_22_TSb <- ape::pcoa(AdNewt_unWUF16Smx_22_TSb)

PCoA_AdNewt_unWUF16Smx_22_TSb_Axes = as.data.frame(PCoA_AdNewt_unWUF16Smx_22_TSb$vectors)
PCoA_AdNewt_unWUF16Smx_22_TSb_Axesb =tibble::rownames_to_column(PCoA_AdNewt_unWUF16Smx_22_TSb_Axes, "SampleID")

PCoA_AdNewt_unWUF16Smx_22_TSb_Axesb %>%
  select(SampleID, Axis.1, Axis.2,Axis.3) %>%
  left_join(AdNewt16S_Meta_Alpha_AntiF, by = "SampleID") %>%
  ggplot(aes(x=Axis.1, y=Axis.2)) +
  geom_point(aes(color=TimeWeek, alpha=0.5,shape =Dose,size =LogCopies_uL))+
  theme_q2r() +
  scale_color_viridis_c(option = "viridis", name="TimeWeek") +
  #stat_ellipse(type = "t",linetype = 2)+
  ggtitle("22C - unweighted")

PCoA_AdNewt_unWUF16Smx_22_TSb_Axesb %>%
  select(SampleID, Axis.1, Axis.2,Axis.3) %>%
  left_join(AdNewt16S_Meta_Alpha_AntiF, by = "SampleID") %>%
  ggplot(aes(x=Axis.1, y=Axis.2,color=Dose,shape =Dose)) +
  geom_point(aes(alpha=0.8,shape =Dose,size =LogCopies_uL))+
  theme_q2r() +
  scale_color_colorblind() +
  #stat_ellipse(type = "t",linetype = 2)+
  ggtitle("22C - unweighted")+
  stat_ellipse(type = "t",linetype = 2) +
  facet_wrap(~TimeWeek)

Final timepoint

AdNewt16S_Meta_Alpha_AntiF_6_TS <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 6) %>%
  filter(Dose == "5x10.3" | Dose == "Control") %>%
  filter(TimeWeek==7) %>%
  drop_na(observed_features)
  
metaList_6_TS<-as.vector(AdNewt16S_Meta_Alpha_AntiF_6_TS$SampleID)

AdNewt_unWUF16Smx2 = AdNewt_unWUF16Smx$data

AdNewt_unWUF16Smx_6_TS <- usedist::dist_subset(AdNewt_unWUF16Smx2,metaList_6_TS)

AdNewt_unWUF16Smx_6_TSb <-as.matrix(AdNewt_unWUF16Smx_6_TS)

AdNewt6_TS_16S_distList<-colnames(AdNewt_unWUF16Smx_6_TSb)

AdNewt16S_Meta_Alpha_AntiF_6_TSb<- AdNewt16S_Meta_Alpha_AntiF_6_TS %>%
  filter(SampleID %in% AdNewt6_TS_16S_distList)

PCoA_AdNewt_unWUF16Smx_6_TSb <- ape::pcoa(AdNewt_unWUF16Smx_6_TSb)

PCoA_AdNewt_unWUF16Smx_6_TSb_Axes = as.data.frame(PCoA_AdNewt_unWUF16Smx_6_TSb$vectors)
PCoA_AdNewt_unWUF16Smx_6_TSb_Axesb =tibble::rownames_to_column(PCoA_AdNewt_unWUF16Smx_6_TSb_Axes, "SampleID")

PCoA_AdNewt_unWUF16Smx_6_TSb_Axesb %>%
  select(SampleID, Axis.1, Axis.2,Axis.3) %>%
  left_join(AdNewt16S_Meta_Alpha_AntiF, by = "SampleID") %>%
  ggplot(aes(x=Axis.1, y=Axis.2,color=Dose)) +
  geom_point(aes(alpha=0.5,shape =Dose,size =LogCopies_uL))+
  theme_q2r() +
  scale_color_viridis_d(option = "viridis", name="TimeWeek") +
  stat_ellipse(type = "t",linetype = 2,level = .95)+
  ggtitle("6C - unweighted")

AdNewt16S_Meta_Alpha_AntiF_14_TS <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 14) %>%
  filter(Dose == "5x10.3" | Dose == "Control") %>%
  filter(TimeWeek==4) %>%
  drop_na(observed_features)
  
metaList_14_TS<-as.vector(AdNewt16S_Meta_Alpha_AntiF_14_TS$SampleID)

AdNewt_unWUF16Smx2 = AdNewt_unWUF16Smx$data

AdNewt_unWUF16Smx_14_TS <- usedist::dist_subset(AdNewt_unWUF16Smx2,metaList_14_TS)

AdNewt_unWUF16Smx_14_TSb <-as.matrix(AdNewt_unWUF16Smx_14_TS)

AdNewt14_TS_16S_distList<-colnames(AdNewt_unWUF16Smx_14_TSb)

AdNewt16S_Meta_Alpha_AntiF_14_TSb<- AdNewt16S_Meta_Alpha_AntiF_14_TS %>%
  filter(SampleID %in% AdNewt14_TS_16S_distList)

PCoA_AdNewt_unWUF16Smx_14_TSb <- ape::pcoa(AdNewt_unWUF16Smx_14_TSb)

PCoA_AdNewt_unWUF16Smx_14_TSb_Axes = as.data.frame(PCoA_AdNewt_unWUF16Smx_14_TSb$vectors)
PCoA_AdNewt_unWUF16Smx_14_TSb_Axesb =tibble::rownames_to_column(PCoA_AdNewt_unWUF16Smx_14_TSb_Axes, "SampleID")

PCoA_AdNewt_unWUF16Smx_14_TSb_Axesb %>%
  select(SampleID, Axis.1, Axis.2,Axis.3) %>%
  left_join(AdNewt16S_Meta_Alpha_AntiF, by = "SampleID") %>%
  ggplot(aes(x=Axis.1, y=Axis.2,color=Dose)) +
  geom_point(aes(alpha=0.5,shape =Dose,size =LogCopies_uL))+
  theme_q2r() +
  scale_color_viridis_d(option = "viridis", name="TimeWeek") +
  stat_ellipse(type = "t",linetype = 2)+
  ggtitle("14C - unweighted")

AdNewt16S_Meta_Alpha_AntiF_22_TS <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 22) %>%
  filter(Dose == "5x10.3" | Dose == "Control") %>%
  filter(TimeWeek==7) %>%
  drop_na(observed_features)
  
metaList_22_TS<-as.vector(AdNewt16S_Meta_Alpha_AntiF_22_TS$SampleID)

AdNewt_unWUF16Smx2 = AdNewt_unWUF16Smx$data

AdNewt_unWUF16Smx_22_TS <- usedist::dist_subset(AdNewt_unWUF16Smx2,metaList_22_TS)

AdNewt_unWUF16Smx_22_TSb <-as.matrix(AdNewt_unWUF16Smx_22_TS)

AdNewt22_TS_16S_distList<-colnames(AdNewt_unWUF16Smx_22_TSb)

AdNewt16S_Meta_Alpha_AntiF_22_TSb<- AdNewt16S_Meta_Alpha_AntiF_22_TS %>%
  filter(SampleID %in% AdNewt22_TS_16S_distList)

PCoA_AdNewt_unWUF16Smx_22_TSb <- ape::pcoa(AdNewt_unWUF16Smx_22_TSb)

PCoA_AdNewt_unWUF16Smx_22_TSb_Axes = as.data.frame(PCoA_AdNewt_unWUF16Smx_22_TSb$vectors)
PCoA_AdNewt_unWUF16Smx_22_TSb_Axesb =tibble::rownames_to_column(PCoA_AdNewt_unWUF16Smx_22_TSb_Axes, "SampleID")

PCoA_AdNewt_unWUF16Smx_22_TSb_Axesb %>%
  select(SampleID, Axis.1, Axis.2,Axis.3) %>%
  left_join(AdNewt16S_Meta_Alpha_AntiF, by = "SampleID") %>%
  ggplot(aes(x=Axis.1, y=Axis.2,color=Dose)) +
  geom_point(aes(alpha=0.5,shape =Dose,size =LogCopies_uL))+
  theme_q2r() +
  scale_color_viridis_d(option = "viridis", name="TimeWeek") +
  stat_ellipse(type = "t",linetype = 2)+
  ggtitle("22C - unweighted")

Taxa

Temp 6 - Week 7

Last week of “good” sample size of both groups

AdNewt16S_Meta_Alpha_AntiF_6_T7 <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 6) %>%
  filter(Dose == "5x10.3" | Dose == "Control") %>%
  filter(TimeWeek==7) %>%
  drop_na(observed_features)
  
metaList_6_T7<-as.vector(AdNewt16S_Meta_Alpha_AntiF_6_T7$SampleID)

AdNewt_UnRarified_Table_T6_w7 = AdNewt_UnRarified_Table %>%
  select(starts_with("NSF6")) %>%
  select(contains("7M"))

AdNewt_UnRarified_Table_T6_w7 = AdNewt_UnRarified_Table_T6_w7[rowSums(AdNewt_UnRarified_Table_T6_w7 != 0)>0,]

otunames <- rownames(AdNewt_UnRarified_Table_T6_w7) 
AdNewt_UnRarified_Table_T6_w7_log = AdNewt_UnRarified_Table_T6_w7 %>% 
  mutate_all(., ~ ifelse(.!=0, log10(.+0.00000000001), log10(1)))

rownames(AdNewt_UnRarified_Table_T6_w7_log)<-otunames

taxasums6w7 = summarize_taxa(AdNewt_UnRarified_Table_T6_w7_log,Taxonomy)$Genus

taxa_heatmap(taxasums6w7, AdNewt16S_Meta_Alpha_AntiF, "Dose",normalize = "none")

Temp 6 - Post exposure across doses

AdNewt16S_Meta_Alpha_AntiF_6_T1 <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 6) %>%
  filter(TimeWeek==1)
  
metaList_6_T1<-as.vector(AdNewt16S_Meta_Alpha_AntiF_6_T1$SampleID)

AdNewt_UnRarified_Table_T6_w1 = AdNewt_UnRarified_Table %>%
  select(starts_with("NSF6")) %>%
  select(contains("-1M"))

AdNewt_UnRarified_Table_T6_w1 = AdNewt_UnRarified_Table_T6_w1[rowSums(AdNewt_UnRarified_Table_T6_w1 != 0)>0,]

 otunames <- rownames(AdNewt_UnRarified_Table_T6_w1) 
AdNewt_UnRarified_Table_T6_w1_log = AdNewt_UnRarified_Table_T6_w1 %>% 
  mutate_all(., ~ ifelse(.!=0, log10(.+0.00000000001), log10(1)))

rownames(AdNewt_UnRarified_Table_T6_w1_log)<-otunames

taxasums = summarize_taxa(AdNewt_UnRarified_Table_T6_w1_log,Taxonomy)$Genus

taxa_heatmap(taxasums, AdNewt16S_Meta_Alpha_AntiF_6_T1, "Dose", normalize = "none")

Temp 14 - Week 4

Last week of “good” sample size of both groups

AdNewt16S_Meta_Alpha_AntiF_14_T4 <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 14) %>%
  filter(Dose == "5x10.3" | Dose == "Control") %>%
  filter(TimeWeek==4) %>%
  drop_na(observed_features)
  
metaList_14_T4<-as.vector(AdNewt16S_Meta_Alpha_AntiF_14_T4$SampleID)

AdNewt_UnRarified_Table_T14_w4 = AdNewt_UnRarified_Table %>%
  select(starts_with("NSF14")) %>%
  select(contains("4M"))

AdNewt_UnRarified_Table_T14_w4 = AdNewt_UnRarified_Table_T14_w4[rowSums(AdNewt_UnRarified_Table_T14_w4 != 0)>0,]

 otunames <- rownames(AdNewt_UnRarified_Table_T14_w4) 
AdNewt_UnRarified_Table_T14_w4_log = AdNewt_UnRarified_Table_T14_w4 %>% 
  mutate_all(., ~ ifelse(.!=0, log10(.+0.00000000001), log10(1)))

rownames(AdNewt_UnRarified_Table_T14_w4_log)<-otunames

taxasums14w4 = summarize_taxa(AdNewt_UnRarified_Table_T14_w4_log,Taxonomy)$Genus

taxa_heatmap(taxasums14w4, AdNewt16S_Meta_Alpha_AntiF, "Dose",normalize = "none")

Temp 14 - Post exposure across doses

AdNewt16S_Meta_Alpha_AntiF_14_T1 <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 14) %>%
  filter(TimeWeek==1) %>%
  drop_na(observed_features)
  
metaList_14_T1<-as.vector(AdNewt16S_Meta_Alpha_AntiF_14_T4$SampleID)

AdNewt_UnRarified_Table_T14_w1 = AdNewt_UnRarified_Table %>%
  select(starts_with("NSF14")) %>%
  select(contains("1M"))

AdNewt_UnRarified_Table_T14_w1 = AdNewt_UnRarified_Table_T14_w1[rowSums(AdNewt_UnRarified_Table_T14_w1 != 0)>0,]

 otunames <- rownames(AdNewt_UnRarified_Table_T14_w1) 
AdNewt_UnRarified_Table_T14_w1_log = AdNewt_UnRarified_Table_T14_w1 %>% 
  mutate_all(., ~ ifelse(.!=0, log10(.+0.00000000001), log10(1)))

rownames(AdNewt_UnRarified_Table_T14_w1_log)<-otunames

taxasums = summarize_taxa(AdNewt_UnRarified_Table_T14_w1_log,Taxonomy)$Genus

taxa_heatmap(taxasums, AdNewt16S_Meta_Alpha_AntiF, "Dose", normalize = "none")

Temp 22 - Week 8

Last week of “good” sample size of both groups

AdNewt16S_Meta_Alpha_AntiF_22_T8 <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 22) %>%
  filter(Dose == "5x10.3" | Dose == "Control") %>%
  filter(TimeWeek==7) %>%
  drop_na(observed_features)
  
metaList_22_T8<-as.vector(AdNewt16S_Meta_Alpha_AntiF_22_T8$SampleID)

AdNewt_UnRarified_Table_T22_w8 = AdNewt_UnRarified_Table %>%
  select(starts_with("NSF22")) %>%
  select(contains("7M"))

AdNewt_UnRarified_Table_T22_w8 = AdNewt_UnRarified_Table_T22_w8[rowSums(AdNewt_UnRarified_Table_T22_w8 != 0)>0,]

 otunames <- rownames(AdNewt_UnRarified_Table_T22_w8) 
AdNewt_UnRarified_Table_T22_w8_log = AdNewt_UnRarified_Table_T22_w8 %>% 
  mutate_all(., ~ ifelse(.!=0, log10(.+0.00000000001), log10(1)))

rownames(AdNewt_UnRarified_Table_T22_w8_log)<-otunames

taxasums22w8 = summarize_taxa(AdNewt_UnRarified_Table_T22_w8_log,Taxonomy)$Genus

taxa_heatmap(taxasums22w8, AdNewt16S_Meta_Alpha_AntiF, "Dose",normalize = "none")

Temp 22 - Post exposure across doses

AdNewt16S_Meta_Alpha_AntiF_22_T1 <- AdNewt16S_Meta_Alpha_AntiF %>%
  filter(Temperature == 22) %>%
  filter(TimeWeek==1) %>%
  drop_na(observed_features)
  
metaList_22_T1<-as.vector(AdNewt16S_Meta_Alpha_AntiF_22_T1$SampleID)

AdNewt_UnRarified_Table_T22_w1 = AdNewt_UnRarified_Table %>%
  select(starts_with("NSF22")) %>%
  select(contains("-1M"))

AdNewt_UnRarified_Table_T22_w1 = AdNewt_UnRarified_Table_T22_w1[rowSums(AdNewt_UnRarified_Table_T22_w1 != 0)>0,]

 otunames <- rownames(AdNewt_UnRarified_Table_T22_w1) 
AdNewt_UnRarified_Table_T22_w1_log = AdNewt_UnRarified_Table_T22_w1 %>% 
  mutate_all(., ~ ifelse(.!=0, log10(.+0.00000000001), log10(1)))

rownames(AdNewt_UnRarified_Table_T22_w1_log)<-otunames

taxasums = summarize_taxa(AdNewt_UnRarified_Table_T22_w1_log,Taxonomy)$Genus

taxa_heatmap(taxasums, AdNewt16S_Meta_Alpha_AntiF, "Dose", normalize = "none")

Differential

QIIME details

### Filtering to Temp 14 and Controls and 10^3 samples

Mollys-MacBook-Pro-2:~ mollybletz$ qiime feature-table filter-samples --i-table /Users/mollybletz/My_FILES/EEID_NSF_Microbiomes/EEID_AqAdult_TempDose_Exp/16S/WithPlasmid/PlasmidWork_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_merged1234/Filtlowsamp_CopCor_True_abund_estimate_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_TempDose_Exp_merged1234.qza --m-metadata-file /Users/mollybletz/My_FILES/EEID_NSF_Microbiomes/EEID_AqAdult_TempDose_Exp/16S/Metadata_NSFEEID_16SRuns_1234Merged_plusExpData.txt --p-where "[Temperature]='14' AND [Dose] IN ('Control','5x10.3')" --o-filtered-table /Users/mollybletz/My_FILES/EEID_NSF_Microbiomes/EEID_AqAdult_TempDose_Exp/16S/WithPlasmid/PlasmidWork_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_merged1234/Temp14-TS_Filtlowsamp_CopCor_True_abund_estimate_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_TempDose_Exp_merged1234.qza --p-no-exclude-ids

### Filtering to OTUs with a cumulative frequency of 100 reads

(qiime2-2020.8) Mollys-MacBook-Pro-2:~ mollybletz$ qiime feature-table filter-features --i-table /Users/mollybletz/My_FILES/EEID_NSF_Microbiomes/EEID_AqAdult_TempDose_Exp/16S/WithPlasmid/PlasmidWork_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_merged1234/Temp14-TS_Filtlowsamp_CopCor_True_abund_estimate_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_TempDose_Exp_merged1234.qza --p-min-frequency 100 --o-filtered-table /Users/mollybletz/My_FILES/EEID_NSF_Microbiomes/EEID_AqAdult_TempDose_Exp/16S/WithPlasmid/PlasmidWork_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_merged1234/Temp14-TS__min100_Filtlowsamp_CopCor_True_abund_estimate_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_TempDose_Exp_merged1234.qza 

###calculating Core features to deetermine threshold

(qiime2-2020.8) Mollys-MacBook-Pro-2:PlasmidWork_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_merged1234 mollybletz$ qiime feature-table core-features --i-table /Users/mollybletz/My_FILES/EEID_NSF_Microbiomes/EEID_AqAdult_TempDose_Exp/16S/WithPlasmid/PlasmidWork_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_merged1234/Temp14-TS__min100_Filtlowsamp_CopCor_True_abund_estimate_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_TempDose_Exp_merged1234.qza --o-visualization Temp14_TS_min100_CoreFeatures.qzv

### Filtering  to OTUs present in at least 20 samples (out of 80 in dataset)

(qiime2-2020.8) Mollys-MacBook-Pro-2:PlasmidWork_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_merged1234 mollybletz$ qiime feature-table filter-features --i-table /Users/mollybletz/My_FILES/EEID_NSF_Microbiomes/EEID_AqAdult_TempDose_Exp/16S/WithPlasmid/PlasmidWork_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_merged1234/Temp14-TS__min100_Filtlowsamp_CopCor_True_abund_estimate_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_TempDose_Exp_merged1234.qza --p-min-samples 20 --o-filtered-table /Users/mollybletz/My_FILES/EEID_NSF_Microbiomes/EEID_AqAdult_TempDose_Exp/16S/WithPlasmid/PlasmidWork_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_merged1234/Temp14-TS_samp20_min100_Filtlowsamp_CopCor_True_abund_estimate_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_TempDose_Exp_merged1234.qza 

Running PCA to identify OTUs that were varying the most across multiple time-points (idea being to reduce the number of OTUS in the mix)

T14_TS_minS20_minr100 = read_qza("./16S/WithPlasmid/PlasmidWork_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_merged1234/Temp14-TS_samp20_min100_Filtlowsamp_CopCor_True_abund_estimate_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_TempDose_Exp_merged1234.qza")

T14_TS_minS20_minr100_phseq = qza_to_phyloseq("./16S/WithPlasmid/PlasmidWork_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_merged1234/Temp14-TS_samp20_min100_Filtlowsamp_CopCor_True_abund_estimate_contamfilt_samplefilt_OTU-table_wplasmid_EEID_AqAdult_TempDose_Exp_merged1234.qza")

T14_TS_minS20_minr100_T1_4 = as.data.frame(T14_TS_minS20_minr100$data) %>%
  select(starts_with("NSF14")) %>%
  select(!contains("5M"))%>%
  select(!contains("6M")) %>%
  select(!contains("7M"))

library(microbiome)


T14_TS_minS20_minr100_T1_4_ps1 = transform(T14_TS_minS20_minr100_T1_4,transform = 'shift', shift = 1)

ord_T14_TS_minS20_minr100_T1_4_ps1_false = prcomp(T14_TS_minS20_minr100_T1_4_ps1, scale. = FALSE)

factoextra::fviz_eig(ord_T14_TS_minS20_minr100_T1_4_ps1_false)

res.var <- factoextra::get_pca_var(ord_T14_TS_minS20_minr100_T1_4_ps1_false)

sample_factorloads = as.data.frame(res.var$contrib )

res.ind <- factoextra::get_pca_ind(ord_T14_TS_minS20_minr100_T1_4_ps1_false)

OTU_factorloads = as.data.frame(res.ind$contrib)

OTU_factorloads_filt2 = OTU_factorloads %>%
  select(Dim.1,Dim.2,Dim.3) %>%
  filter(Dim.1>0.1 | Dim.2>0.1| Dim.3>0.1)

PCA_OTUS = rownames(OTU_factorloads_filt2)

Ploting OTUs through time

T14_TS_minS20_minr100_1_4_PCAotus = T14_TS_minS20_minr100_T1_4_ps1 %>%
  rownames_to_column(., var = "OTUID") %>%
  filter(OTUID %in% PCA_OTUS)%>%
  column_to_rownames(., var = "OTUID")

AdNewt16S_Meta_Alpha_AntiF_time_load = AdNewt16S_Meta_Alpha_AntiF %>%
 select(SampleID,TimeWeek,LogCopies_uL,Dose,IndividualID)

T14_TS_minS20_minr100_1_4_PCAotus_tr_meta = t(T14_TS_minS20_minr100_1_4_PCAotus) %>%
  as.data.frame() %>%
  rownames_to_column(., var = "SampleID") %>%
  pivot_longer(cols = 2:68, names_to = "OTUID", values_to = "abundance") %>%
  left_join(.,AdNewt16S_Meta_Alpha_AntiF_time_load, by = "SampleID")

  T14_TS_minS20_minr100_1_4_PCAotus_tr_meta$IndividualID = as.factor(T14_TS_minS20_minr100_1_4_PCAotus_tr_meta$IndividualID)

  
T14_TS_minS20_minr100_1_4_PCAotus_tr_meta %>%
  ggplot(aes(x = TimeWeek, y = log10(as.numeric(abundance+1)),color = Dose ))+
  stat_summary(geom="errorbar", fun.data=mean_se, width=0) +
  stat_summary(geom="point", fun.data=mean_se) +
  stat_summary(geom="line", fun.data=mean_se) +
  xlab("Week") +
  ylab(" OTU abundance") +
  theme_q2r() + 
  scale_color_manual(values = c("#D37105","#475C63")) +
  facet_wrap(~OTUID, scales = "free")

modeling explorations

#models

T14_TS_minS20_minr100_1_4_PCAotus_tr_meta_OTU1 = T14_TS_minS20_minr100_1_4_PCAotus_tr_meta %>%
  filter(OTUID== "affbf2c1be9ca467eb201941d9782b47")

hist(T14_TS_minS20_minr100_1_4_PCAotus_tr_meta_OTU1$abundance)

mod = glmer(ceiling(abundance) ~ TimeWeek + Dose + TimeWeek*Dose + (1|IndividualID), family = poisson,data = T14_TS_minS20_minr100_1_4_PCAotus_tr_meta_OTU1)

summary(mod)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: ceiling(abundance) ~ TimeWeek + Dose + TimeWeek * Dose + (1 |  
##     IndividualID)
##    Data: T14_TS_minS20_minr100_1_4_PCAotus_tr_meta_OTU1
## 
##       AIC       BIC    logLik  deviance  df.resid 
##  803582.8  803593.7 -401786.4  803572.8        60 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -178.47  -44.66  -21.62   -0.03  543.75 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  IndividualID (Intercept) 5.923    2.434   
## Number of obs: 65, groups:  IndividualID, 15
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         5.698736   1.033208   5.516 3.48e-08 ***
## TimeWeek            0.350731   0.002343 149.690  < 2e-16 ***
## Dose5x10.3          0.938237   1.268220   0.740    0.459    
## TimeWeek:Dose5x10.3 0.005062   0.002627   1.927    0.054 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TimeWk D510.3
## TimeWeek    -0.006              
## Dose5x10.3  -0.800  0.005       
## TmWk:D510.3  0.005 -0.892 -0.005
## convergence code: 0
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?

T14_TS_minS20_minr100_1_4_PCAotus_tr_meta_OTU2 = T14_TS_minS20_minr100_1_4_PCAotus_tr_meta %>%
  filter(OTUID== "070cc7cefe482682f056e340a569c569")

hist(T14_TS_minS20_minr100_1_4_PCAotus_tr_meta_OTU2$abundance)

mod1 = glmer(ceiling(abundance) ~ TimeWeek + Dose + TimeWeek*Dose + (1|IndividualID), family = poisson,data = T14_TS_minS20_minr100_1_4_PCAotus_tr_meta_OTU2)

summary(mod1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: ceiling(abundance) ~ TimeWeek + Dose + TimeWeek * Dose + (1 |  
##     IndividualID)
##    Data: T14_TS_minS20_minr100_1_4_PCAotus_tr_meta_OTU2
## 
##      AIC      BIC   logLik deviance df.resid 
## 114033.7 114044.6 -57011.9 114023.7       60 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -137.165   -3.126   -0.552    0.469  214.336 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  IndividualID (Intercept) 6.302    2.51    
## Number of obs: 65, groups:  IndividualID, 15
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          0.50649    1.14347   0.443   0.6578    
## TimeWeek             0.22683    0.05620   4.036 5.44e-05 ***
## Dose5x10.3           2.46395    1.39398   1.768   0.0771 .  
## TimeWeek:Dose5x10.3  0.02575    0.05628   0.458   0.6473    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TimeWk D510.3
## TimeWeek    -0.129              
## Dose5x10.3  -0.820  0.106       
## TmWk:D510.3  0.129 -0.999 -0.106

## Group model

glmermods = groupedstats::grouped_glmer(T14_TS_minS20_minr100_1_4_PCAotus_tr_meta, grouping.vars = "OTUID", floor(abundance) ~ TimeWeek + Dose + TimeWeek*Dose + (1|IndividualID), family = poisson, output = "tidy",tidy.args = list(conf.int = TRUE, conf.level = 0.95, effects = "fixed", conf.method ="Wald")) 


glmermods %>%
  ggplot(aes(x = term, y = statistic)) +
  geom_point(aes(color = p.value < 0.01))+
  coord_flip() +
  facet_wrap(~OTUID)

glmermods_wide = glmermods %>%
  select(OTUID,p.value,term) %>%
  pivot_wider(id_cols = "OTUID",names_from = term, values_from  = p.value)

glmermods_wide_sig_Dose = glmermods_wide %>%
  filter(Dose5x10.3<0.001) %>%
  filter(!`TimeWeek:Dose5x10.3`<0.001)%>%
  filter(!TimeWeek<0.001)

glmermods_wide_sig_Inter = glmermods_wide %>%
  filter(`TimeWeek:Dose5x10.3`<0.001)

SigOTUs_T14 = glmermods_wide_sig_Inter$OTUID
  

SigOTUs_T14 = as.data.frame(SigOTUs_T14) %>%
  add_column(Sig="Yes") %>%
  rename(OTUID=SigOTUs_T14)
  
  
T14_TS_minS20_minr100_1_4_PCAotus_tr_meta %>%
  left_join(SigOTUs_T14,by =  "OTUID") %>%
  mutate(Sig = replace_na(Sig,"No")) %>%
  ggplot(aes(x = TimeWeek, y = log10(as.numeric(abundance+1)),color = Dose))+
  stat_summary(geom="errorbar", fun.data=mean_se, width=0) +
  stat_summary(geom="point", fun.data=mean_se) +
  stat_summary(geom="line", fun.data=mean_se) +
  xlab("Week") +
  ylab(" log OTU abundance") +
  theme_q2r() + 
  scale_color_manual(values = c("#D37105","#475C63")) +
  facet_wrap(~OTUID, scales = "free")

ALL otus not PCA filtered

T14_TS_minS20_minr100_1_4_tr_meta = T14_TS_minS20_minr100_T1_4 %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column(., var = "SampleID") %>%
  pivot_longer(cols = 2:72, names_to = "OTUID", values_to = "abundance") %>%
  left_join(.,AdNewt16S_Meta_Alpha_AntiF_time_load, by = "SampleID")

T14_TS_minS20_minr100_1_4_tr_meta %>%
  ggplot(aes(x = TimeWeek, y = log10(as.numeric(abundance+1)),color = Dose))+
  stat_summary(geom="errorbar", fun.data=mean_se, width=0) +
  stat_summary(geom="point", fun.data=mean_se) +
  stat_summary(geom="line", fun.data=mean_se) +
  xlab("Week") +
  ylab(" log OTU abundance") +
  theme_q2r() + 
  scale_color_manual(values = c("#D37105","#475C63")) +
  facet_wrap(~OTUID, scales = "free")

otunames <- rownames(T14_TS_minS20_minr100_1_4_PCAotus) 

T14_TS_minS20_minr100_1_4_PCAotus_log = T14_TS_minS20_minr100_1_4_PCAotus %>% 
  mutate_all(., ~ ifelse(.!=0, log10(.+0.00000000001), log10(1)))
rownames(T14_TS_minS20_minr100_1_4_PCAotus_log)<-otunames

taxasums14TS = summarize_taxa(T14_TS_minS20_minr100_1_4_PCAotus_log,Taxonomy)$Genus

taxa_heatmap(taxasums14TS, AdNewt16S_Meta_Alpha_AntiF, "DoseTime",normalize = "none")